## Greener Governments: Partisan Ideologies, Executive Institutions and Environmental Policies." Forthcoming in Environmental Politics. (with R. Thomson)
# Replication file: 2016-01-02


library(foreign)
library(car)
library(stargazer)
library(nnet)
library(simex)
library(lme4)


# Import data set (includes data from CMP, DEU, RAI and controls specified in the paper
#d.all <- read.csv("All-Model-Relevant-Data.csv") 

# 1,093 obs on all created vars, restrict to 1,022?
# Remove unused vars

#d.all2 <- d.all[,c("envpos.ord", "logCMP.PM", "logCMP.ENV", "loggdppc", "countp2ms", "com.100.env", "ep.100.env", "Kasscent", "envpos", "msid", "prnrnmc", "isnrnmc", "abbrev", "COW.Code", "inyr", "envpos.0.1", "envpos.0.2", "CMP.ENV.ENV", "CMP.ENV.PM", "CMP.ENV.PM.SE", "CMP.ENV.ENV.SE", "logCMP.PM.SE", "logCMP.ENV.SE")]

#write.csv(d.all2, "All-Model-Relevant-Data2.csv") 

d.all <- read.csv("All-Model-Relevant-Data2.csv") 


#####################################################
######## Regressions
all.1 <- multinom(envpos.ord~logCMP.PM+logCMP.ENV+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.all)
dim(model.frame(all.1))

all.2 <- multinom(envpos.ord~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.all)
dim(model.frame(all.2))

d.same <- subset(d.all, d.all$logCMP.PM==d.all$logCMP.ENV)
same.3 <- multinom(envpos.ord~logCMP.PM+Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.same)
dim(model.frame(same.3))

d.different <- subset(d.all, d.all$logCMP.PM!=d.all$logCMP.ENV)
different.4 <- multinom(envpos.ord~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.different)
dim(model.frame(different.4))


############################################################
# TABLE 2
stargazer(all.1, all.2, same.3, different.4, digits=2, star.cutoffs=c(.05,.01,.001))


#################################
# FIG 1: Coefficient dot plots ("The Importance of Party Ideology and Institutional Context for Environmental Position-taking")
# Bootstrap CIs
# NOTE: Random sampling will draw different CIs each time

# Focus on model observations only
d.pop <- model.frame(all.2) # 909 obs

# Create random samples
reps <- 1000
sampling.obj <- matrix(nrow = dim(d.pop)[1], ncol = reps, data = NA)

for (i in 1:reps){
  sampling.obj[,i] <- sample(1:dim(d.pop)[1], dim(d.pop)[1], replace=T)
}

# Store coefs for 0 vs 1 and 0 vs 2 for both actors, just related to marginal effects
par.est <- matrix(nrow=reps, ncol=8, data=NA)

for (i in 1:reps){
  # run multinom on d.pop ordered by one column of the sampling object
  reg.loop <- multinom(envpos.ord~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.pop[sampling.obj[,i],])
  # Save coefficients
  par.est[i,] <- coef(reg.loop)[c(3,17,7,19,4,18,8,20)]
}

head(par.est)
colnames(par.est) <- c("PM.1", "PM.x.1", "ENV.1", "ENV.x.1", "PM.2", "PM.x.2", "ENV.2", "ENV.x.2")
head(par.est)

apply(X=par.est,  MARGIN=2, FUN=quantile, probs=.05)
apply(X=par.est,  MARGIN=2, FUN=mean)
apply(X=par.est,  MARGIN=2, FUN=quantile, probs=.95)

# 0 vs 1, PM central, decentral, diff, ENV central, decentral, diff 
quantile(par.est[,1]+par.est[,2], probs = c(.05,.5,.95))
quantile(par.est[,1], probs = c(.05,.5,.95))
quantile(par.est[,2], probs = c(.05,.5,.95))

quantile(par.est[,3]+par.est[,4], probs = c(.05,.5,.95))
quantile(par.est[,3], probs = c(.05,.5,.95))
quantile(par.est[,4], probs = c(.05,.5,.95))


# 0 vs 2, PM central, decentral, diff, ENV central, decentral, diff 
quantile(par.est[,5]+par.est[,6], probs = c(.05,.5,.95))
quantile(par.est[,5], probs = c(.05,.5,.95))
quantile(par.est[,6], probs = c(.05,.5,.95))

quantile(par.est[,7]+par.est[,8], probs = c(.05,.5,.95))
quantile(par.est[,7], probs = c(.05,.5,.95))
quantile(par.est[,8], probs = c(.05,.5,.95))


# Save effects and change differences into absolute value! Or don't?
points.1 <- c(quantile(par.est[,1]+par.est[,2], probs = .5), quantile(par.est[,1], probs=.5), quantile(par.est[,2], probs =.5), quantile(par.est[,3]+par.est[,4], probs =.5), quantile(par.est[,3], probs =.5), quantile(par.est[,4], probs =.5)*-1)

low.1 <- c(quantile(par.est[,1]+par.est[,2], probs = .025), quantile(par.est[,1], probs=.025), quantile(par.est[,2], probs =.025), quantile(par.est[,3]+par.est[,4], probs =.025), quantile(par.est[,3], probs =.025), quantile(par.est[,4], probs =.025)*-1)

high.1 <- c(quantile(par.est[,1]+par.est[,2], probs = .975), quantile(par.est[,1], probs=.975), quantile(par.est[,2], probs =.975), quantile(par.est[,3]+par.est[,4], probs =.975), quantile(par.est[,3], probs =.975), quantile(par.est[,4], probs =.975)*-1)

groups <- c(rep("PM",3),rep("Env Minister",3))
tags <- rep(c("Centralized","Decentralized","Abs Difference"),2)

f.1 <- data.frame(groups, tags, low.1, points.1, high.1)
f.1

levels(f.1$groups)
f.1$groups <- relevel(f.1$groups, "PM")
f.1

# 0 vs 1
dotchart(f.1$points.1[6:1], groups=f.1$groups[6:1], labels=f.1$tags[6:1], xlim=c(-1,2), pch=19)
abline(v=0, lty=2, col="grey")
segments(f.1$low.1, c(8,7,6,3,2,1), f.1$high.1, c(8,7,6,3,2,1))


# 0 vs 2
points <- c(quantile(par.est[,5]+par.est[,6], probs = .5), quantile(par.est[,5], probs=.5), quantile(par.est[,6], probs =.5), quantile(par.est[,7]+par.est[,8], probs =.5), quantile(par.est[,7], probs =.5), quantile(par.est[,8], probs =.5)*-1)

low <- c(quantile(par.est[,5]+par.est[,6], probs = .025), quantile(par.est[,5], probs=.025), quantile(par.est[,6], probs =.025), quantile(par.est[,7]+par.est[,8], probs =.025), quantile(par.est[,7], probs =.025), quantile(par.est[,8], probs =.025)*-1)

high <- c(quantile(par.est[,5]+par.est[,6], probs = .975), quantile(par.est[,5], probs=.975), quantile(par.est[,6], probs =.975), quantile(par.est[,7]+par.est[,8], probs =.975), quantile(par.est[,7], probs =.975), quantile(par.est[,8], probs =.975)*-1)

groups <- c(rep("PM",3),rep("Env Minister",3))
tags <- rep(c("Centralized","Decentralized","Abs Difference"),2)

f <- data.frame(groups, tags, low, points, high)
f

levels(f$groups)
f$groups <- relevel(f$groups, "PM")
f

#png("../../Writing/CURRENT_VERSION/Figures_Tables/Coef_Effects.png")
dotchart(f$points[6:1], groups=f$groups[6:1], labels=f$tags[6:1], xlim=c(-1,2), pch=19)
abline(v=0, lty=2, col="grey")
segments(f$low, c(8,7,6,3,2,1), f$high, c(8,7,6,3,2,1))
#dev.off()



##########################################################
### FIG 2: Predicted probs - envenv + decentral, all else at mean
# Bootstrap CIs

# Restrict data to only complete rows
d.pop <- model.frame(all.2) # 909 obs

# Sample from d.pop with replacement
sample(1:909, 909, replace=T)

# Repeat sampling 1,000 times, save all in object as different columns
sampling.obj <- matrix(nrow = 909, ncol = 1000,data = NA)

for (i in 1:1000){
  sampling.obj[,i] <- sample(1:909, 909, replace=T)
}

newdat.PM <- data.frame(logCMP.PM=seq(0,2.6, by=.02), logCMP.ENV=rep(mean(d.all$logCMP.ENV, na.rm=T), 131), Kasscent=rep(1,131), loggdppc=rep(mean(d.all$loggdppc, na.rm=T), 131), countp2ms=rep(mean(d.all$countp2ms, na.rm=T), 131), com.100.env=rep(0, 131), ep.100.env=rep(1, 131))

newdat.ENV <- data.frame(logCMP.PM=rep(mean(d.all$logCMP.PM, na.rm=T), 77), logCMP.ENV=seq(0,3.812, by=.05), Kasscent=rep(0,77), loggdppc=rep(median(d.all$loggdppc, na.rm=T), 77), countp2ms=rep(mean(d.all$countp2ms, na.rm=T), 77), com.100.env=rep(0, 77), ep.100.env=rep(1, 77))

# Create objects to store predicted probs
PM.Strong <- matrix(nrow=131, ncol = 1000, data = NA)
PM.Weak <- matrix(nrow=131, ncol = 1000, data = NA)
ENV.Strong <- matrix(nrow=77, ncol = 1000, data = NA)
ENV.Weak <- matrix(nrow=77, ncol = 1000, data = NA)


for (i in 1:1000){
  # run multinom on d.pop ordered by one column of the sampling object
  reg.loop <- multinom(envpos.ord~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.pop[sampling.obj[,i],])
  
  # Save predicted probs for PM in central - weakest and strongest
  newdat.pred <- predict(reg.loop, newdat.PM, type ="probs")
  PM.Weak[,i] <- newdat.pred[,1]
  PM.Strong[,i] <- newdat.pred[,3]
  
  # Save predicted probs ENV in decentral - weakest and strongest
  newdat.pred <- predict(reg.loop, newdat.ENV, type ="probs")
  ENV.Weak[,i] <- newdat.pred[,1]
  ENV.Strong[,i] <- newdat.pred[,3]
}


## Make plots with simulated probability predictions
# Mean of each row as solid line
# 90% CI: 5% prediction line and 95% prediction line

apply(X = ENV.Weak, MARGIN = 1, FUN = mean)
apply(X = ENV.Weak, MARGIN = 1, quantile, probs=.05)
apply(X = ENV.Weak, MARGIN = 1, quantile, probs=.95)

# ENV MINISTER

#png("../../Writing/CURRENT_VERSION/Figures_Tables/EnvMinister-Pred-Prob-CMP-CIs.png")
a <- c(1,2,4,8,16,32)
plot(seq(0,3.812, by=.05), apply(X = ENV.Weak, MARGIN = 1, FUN = mean), lwd=4, ylim=c(0,1), xlim=c(0,4), type="l", xlab="log(CMP Environment Minister + 1)", ylab="Probability", yaxt="n", cex.axis=1.3, cex.lab=1.3, col="grey")
axis(2, at=seq(0,1,.2), labels=c("0.0", "0.2", "0.4", "0.6", "0.8", "1.0"), las=1, tick=F, cex.axis=1.3)
par(new=T)
plot(seq(0,3.812, by=.05), apply(X = ENV.Weak, MARGIN = 1, quantile, probs=.05), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="", col="grey")
par(new=T)
plot(seq(0,3.812, by=.05), apply(X = ENV.Weak, MARGIN = 1, quantile, probs=.95), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="", col="grey")
text(0, .87, "Weakest", pos=4, cex=1.2, col="dimgrey")
axis(3, at=log(a+1), labels=a, cex.axis=1.3, lwd=0, lwd.ticks = 1)
mtext(text = "CMP Environment Minister (%)", side = 3, line = 2.5, cex = 1.3)
par(new=T)
plot(seq(0,3.812, by=.05), apply(X = ENV.Strong, MARGIN = 1, FUN = mean), lwd=4, ylim=c(0,1), xlim=c(0,4), type="l", xlab="", ylab="", yaxt="n", cex.axis=1.3, cex.lab=1.3)
par(new=T)
plot(seq(0,3.812, by=.05), apply(X = ENV.Strong, MARGIN = 1, quantile, probs=.05), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="")
par(new=T)
plot(seq(0,3.812, by=.05), apply(X = ENV.Strong, MARGIN = 1, quantile, probs=.95), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="")
text(0, .17, "Strongest", pos=4, cex=1.2)
abline(h=1, lty=2, col="gray")
abline(h=.8, lty=2, col="gray")
abline(h=.6, lty=2, col="gray")
abline(h=.4, lty=2, col="gray")
abline(h=.2, lty=2, col="gray")
#dev.off()


# PM

#png("../../Writing/CURRENT_VERSION/Figures_Tables/PM-Pred-Prob-CMP-CIs.png")
a <- c(1,2,4,8,16,32)
plot(seq(0,2.6, by=.02), apply(X = PM.Weak, MARGIN = 1, FUN = mean), lwd=4, ylim=c(0,1), xlim=c(0,4), type="l", xlab="log(CMP PM + 1)", ylab="Probability", yaxt="n", cex.axis=1.3, cex.lab=1.3, col="grey")
axis(2, at=seq(0,1,.2), labels=c("0.0", "0.2", "0.4", "0.6", "0.8", "1.0"), las=1, tick=F, cex.axis=1.3, col="grey")
par(new=T)
plot(seq(0,2.6, by=.02), apply(X = PM.Weak, MARGIN = 1, quantile, probs=.05), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="", col="grey")
par(new=T)
plot(seq(0,2.6, by=.02), apply(X = PM.Weak, MARGIN = 1, quantile, probs=.95), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="", col="grey")
text(0, .78, "Weakest", pos=4, cex=1.2, col="dimgrey")
axis(3, at=log(a+1), labels=a, cex.axis=1.3, lwd=0, lwd.ticks = 1)
mtext(text = "CMP PM (%)", side = 3, line = 2.5, cex = 1.3)
par(new=T)
plot(seq(0,2.6, by=.02), apply(X = PM.Strong, MARGIN = 1, FUN = mean), lwd=4, ylim=c(0,1), xlim=c(0,4), type="l", xlab="", ylab="", yaxt="n", cex.axis=1.3, cex.lab=1.3)
par(new=T)
plot(seq(0,2.6, by=.02), apply(X = PM.Strong, MARGIN = 1, quantile, probs=.05), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="")
par(new=T)
plot(seq(0,2.6, by=.02), apply(X = PM.Strong, MARGIN = 1, quantile, probs=.95), lwd=1, lty=2, ylim=c(0,1), xlim=c(0,4), type="l", bty="n", yaxt="n", xaxt="n", ylab="", xlab="")
text(0, .23, "Strongest", pos=4, cex=1.2)
abline(h=1, lty=2, col="gray")
abline(h=.8, lty=2, col="gray")
abline(h=.6, lty=2, col="gray")
abline(h=.4, lty=2, col="gray")
abline(h=.2, lty=2, col="gray")
#dev.off()



##############
## APPENDIX ##


###################################
## Table A.2: Summary Statistics
d.summary <- model.frame(all.2)

# To get predictors not included in main regression (e.g. envpos, CMP.ENV.ENV, CMP.ENV.PM), match cases using those vars
d.test <- d.all[complete.cases(d.all$envpos, d.all$CMP.ENV.PM, d.all$Kasscent, d.all$CMP.ENV.ENV, d.all$loggdppc, d.all$countp2ms, d.all$com.100.env, d.all$ep.100.env),]
dim(d.test)

summary(d.test$envpos)
sd(d.test$envpos)
table(d.summary$envpos.ord)
table(d.summary$envpos.ord)/length(d.summary$envpos.ord)

summary(d.test$CMP.ENV.PM)
sd(d.test$CMP.ENV.PM)
summary(d.summary$logCMP.PM)
sd(d.summary$logCMP.PM)

summary(d.test$CMP.ENV.ENV)
sd(d.test$CMP.ENV.ENV)
summary(d.summary$logCMP.ENV)
sd(d.summary$logCMP.ENV)

summary(d.summary$Kasscent)
sd(d.summary$Kasscent)
summary(d.summary$loggdppc)
sd(d.summary$loggdppc)
summary(d.summary$countp2ms)
sd(d.summary$countp2ms)
mean(d.summary$com.100.env)
mean(d.summary$ep.100.env)


###############################
## FIG A.1: DV construction
# a. Histogram of d.narrow$envpos
#jpeg("../Figures_Tables/Envpos-Hist.jpg")
hist(d.all$envpos[!is.na(d.all$msid)], main="", col="gray", ylim=c(0,500), cex=1.5, cex.axis=1.5, cex.lab=1.5, xlab="Environmental Position", ylab="", las=1)
#dev.off()

# b. Barplot of d.narrow$envpos.ord
#jpeg("../Figures_Tables/Envpos-ord-bars.jpg")
barplot(table(d.all$envpos.ord[!is.na(d.all$msid)]), main="", col="gray", ylim=c(0,500), cex=1.3, cex.axis=1.5, cex.lab=1.5, space=c(.85,.85,.85), xlab="Environmental Position", names.arg=c("Weakest", "Intermediate", "Strongest"), xlim=c(0,7), las=1)
#dev.off()



#################################
## Table A.3: Issue Heterogeneity Robustness Check

# 3. Random effect for issue, Repeat multilevel separately to account for 53% of cases env=pm

res0.1 <- glmer(envpos.0.1~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env+(1|isnrnmc), data=d.all, family=binomial(link="logit"), control=glmerControl(optimizer="bobyqa")) #
#ranef(res0.1)
summary(res0.1)

res0.2 <- glmer(envpos.0.2~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env+(1|isnrnmc), data=d.all, family=binomial(link="logit"), control=glmerControl(optimizer="bobyqa")) #
#ranef(res0.2)
summary(res0.2)

stargazer(res0.1, res0.2, digits=2, star.cutoffs=c(.10,.05,.01), star.char=c(".","*","**"))



#################################
## Table A.4: CMP Measurement Error Robustness Check
# NOTE: Results will differ slightly each simulation

# 0 vs 1
model1.logit.1 <- glm(envpos.0.1~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.all, family=binomial, x=T)
summary(model1.logit.1)

refined.sample.all.2 <- complete.cases(d.all$envpos.0.1, d.all$logCMP.PM, d.all$logCMP.ENV, d.all$Kasscent, d.all$loggdppc, d.all$countp2ms, d.all$com.100.env, d.all$ep.100.env)
table(refined.sample.all.2) # 717 cases

a1 <- simex(model=model1.logit.1, SIMEXvariable=cbind("logCMP.PM","logCMP.ENV"), measurement.error=cbind(d.all$logCMP.PM.SE[refined.sample.all.2], d.all$logCMP.ENV.SE[refined.sample.all.2]), B=1000, asymptotic=F)
summary(a1)


# 0 vs 2
model1.logit.2 <- glm(envpos.0.2~logCMP.PM*Kasscent+logCMP.ENV*Kasscent+loggdppc+countp2ms+com.100.env+ep.100.env, data=d.all, family=binomial, x=T)
summary(model1.logit.2)

refined.sample.all.2 <- complete.cases(d.all$envpos.0.2, d.all$logCMP.PM, d.all$logCMP.ENV, d.all$Kasscent, d.all$loggdppc, d.all$countp2ms, d.all$com.100.env, d.all$ep.100.env)
table(refined.sample.all.2) # 563 cases

a2 <- simex(model=model1.logit.2, SIMEXvariable=cbind("logCMP.PM","logCMP.ENV"), measurement.error=cbind(d.all$logCMP.PM.SE[refined.sample.all.2], d.all$logCMP.ENV.SE[refined.sample.all.2]), B=1000, asymptotic=F)
summary(a2)

summary(a1)
logLik(model1.logit.1)
summary(a2)
logLik(model1.logit.2)


#############################################################
# ROBUSTNESS: Check main model for robustness of Kassim's centralization by repeatedly RECLASSIFYING states and re-running regression.

d.error <- d.all
table(d.error$msid, d.error$Kasscent) # Every state has only one Kasscent score for this range of years

coef.robustness <- matrix(NA, 29, 4) # Coef results object
p.robustness <- matrix(NA, 29, 4) # p values object

for (i in 3:29){
  d.error <- d.all
  d.error$Kasscent2 <- ifelse(d.error$msid==i & d.error$Kasscent==0, 1, ifelse(d.error$msid==i & d.error$Kasscent==1,0, d.error$Kasscent))
  dec.multi <- multinom(envpos.ord~logCMP.PM*Kasscent2+logCMP.ENV*Kasscent2, data=d.error)
  z <- coef(dec.multi)[2,c(2,5,4,6)]/summary(dec.multi)$standard.errors[2,c(2,5,4,6)]
  p.robustness[i,] <- (1-pnorm(abs(z), 0, 1))*2 # 2-tailed z test
  coef.robustness[i,] <- coef(dec.multi)[2,c(2,5,4,6)]
}

coef.robustness
p.robustness

table(p.robustness[,1]<.05) # logCMP.PM 26 insignificant, 1 significant
table(p.robustness[,2]<.05) # xCentral 24 significant, 3 insignificant
table(p.robustness[,3]<.05) # logCMP.ENV All significant
table(p.robustness[,4]<.05) # xCentral 26 significant, 1 insignificant


#############################################################
# ROBUSTNESS: Check model all.0 for robustness by repeatedly dropping states and re-running regression.

all.0 <- multinom(envpos.ord~logCMP.PM*Kasscent+logCMP.ENV*Kasscent, data=d.all)

summary(all.0)
names(all.0)
str(summary(all.0))
coef(all.0)
coef(all.0)[2,c(2,5,4,6)] # CMP.PM, xKasscent, CMP.ENV, xKasscent
summary(all.0)$standard.errors[2,c(2,5,4,6)] # SEs

z <- coef(all.0)[2,c(2,5,4,6)]/summary(all.0)$standard.errors[2,c(2,5,4,6)]
p <- (1-pnorm(abs(z), 0, 1))*2 # 2-tailed z test
p

rbind(coef(all.0)[2,c(2,5,4,6)], p) #

coef.robustness <- matrix(NA, 29, 4) # Coef results object
p.robustness <- matrix(NA, 29, 4) # p values object

for (i in 3:29){
  dec.multi <- multinom(envpos.ord~logCMP.PM*Kasscent+logCMP.ENV*Kasscent, data=d.all, subset=d.all$msid!=i)
  z <- coef(dec.multi)[2,c(2,5,4,6)]/summary(dec.multi)$standard.errors[2,c(2,5,4,6)]
  p.robustness[i,] <- (1-pnorm(abs(z), 0, 1))*2 # 2-tailed z test
  coef.robustness[i,] <- coef(dec.multi)[2,c(2,5,4,6)]
}

coef.robustness
p.robustness

table(p.robustness[,1]<.05) # All insig
table(p.robustness[,2]<.05) # 26 significant, 1 insignificant
table(p.robustness[,3]<.05) # All significant
table(p.robustness[,4]<.05) # All significant